###################################################
### Load survival package and read data
###################################################

library(survival)
eg1a <- read.table("eg1a.dat",header=TRUE)
eg1b <- read.table("eg1b.dat",header=TRUE)


###################################################
### Fit a single Kaplan-Meier curve
###################################################

eg1a.km <- survfit(Surv(days,status),data=eg1a)


###################################################
### Plot a Kaplan-Meier curve
###################################################

plot(eg1a.km,mark.time=FALSE,xlab="Days",ylab="Proportion Event-Free",main="Kaplan-Meier Survivor Function\nwith Greenwood 95% Confidence Interval")


###################################################
### Fit multiple (two here) KM curves on the basis
### of a covariate (rx=treatment here)
###################################################

eg1b.km <- survfit(Surv(days,status)~rx,data=eg1b)


###################################################
### Plot KM for multiple curves
###################################################

plot(eg1b.km,mark.time=FALSE,lty=c(1,2),xlab="Days",ylab="Proportion Event-Free")


###################################################
### Compare survival
###################################################

survdiff(Surv(days,status)~rx,data=eg1b)


###################################################
### Direct calculation of rho for exponential
###################################################

rhohat <- sum(eg1a$status)/sum(eg1a$days)
se.rhohat <- sqrt(rhohat^2 / sum(eg1a$status))
symci <- rhohat + c(-1,1)*qnorm(0.975)*se.rhohat
c(rhohat,symci)


###################################################
### Plot of H(t) = -log(F(t)) vs. t
###################################################

plot(survfit(Surv(days,status)~1,data=eg1a),mark.time=FALSE,conf.int=FALSE,fun="cumhaz",ylab="H(t)",xlab="Time",main="Cumulative Hazard Plot")


###################################################
### Computation of likelihood based CI in exponential
###################################################

llik <- function(rho,events,riskt) {
  events * log(rho) - rho * riskt
}

lrtf <- function(x,rho,events,riskt,alpha=0.05) {
  2*(llik(rho,events,riskt) - llik(x,events,riskt)) - qchisq(1-alpha,1)
}

lrtci <- c(uniroot(lrtf,interval=c(rhohat-3*se.rhohat,rhohat),rho=rhohat,
                   events=sum(eg1a$status),riskt=sum(eg1a$days))$root,
           uniroot(lrtf,interval=c(rhohat,rhohat+3*se.rhohat),rho=rhohat,
                   events=sum(eg1a$status),riskt=sum(eg1a$days))$root)
lrtci


###################################################
### Plot of the exponential log-likelihood
###################################################

rhos <- c(seq(rhohat-3*se.rhohat,rhohat,length=100),seq(rhohat,rhohat+4*se.rhohat,length=100))
pllik <- sapply(rhos,llik,events=sum(eg1a$status),riskt=sum(eg1a$days))

plot(pllik~rhos,type='l',xlab=expression(rho),ylab="log likelihood",main="Exponential Log Likelihood")


###################################################
### Estimate rho by exponential parametric survival model
###################################################

eg1a.sr <- survreg(Surv(days,status)~1,data=eg1a,dist="exponential")
summary(eg1a.sr)
exp(-coef(eg1a.sr))
regrci <- c(exp(-1*(coef(eg1a.sr)+qnorm(0.975)*sqrt(vcov(eg1a.sr))[1,1])),
            exp(-1*(coef(eg1a.sr)-qnorm(0.975)*sqrt(vcov(eg1a.sr))[1,1])))
regrci



###################################################
### Direct hazard ratio estimation for exponential
###################################################

rhohats <- tapply(eg1b$status,eg1b$rx,sum)/tapply(eg1b$days,eg1b$rx,sum)
hr <- rhohats[2]/rhohats[1]
lhrse <- sqrt(sum(1/tapply(eg1b$status,eg1b$rx,sum)))
c(hr, log(hr), lhrse)


###################################################
### Model-based estimation of hazard ratio
###################################################

eg1b.sr <- survreg(Surv(days,status)~rx,data=eg1b,dist="exponential")
summary(eg1b.sr)
exp(-coef(eg1b.sr)[2])


###################################################
### Plot of log(-log(F(t))) vs. log(t) to check for weibull
###################################################

plot(survfit(Surv(days,status)~rx,data=eg1b),mark.time=FALSE,col=c(1,2),
     fun="cloglog",xlim=range(eg1b$days)*c(1,1.5),ylab="log(H(t)",xlab="log(t)")


###################################################
### Fit weibull model
###################################################

eg2b.sr <- survreg(Surv(days,status)~rx,data=eg1b,dist="weibull")
summary(eg2b.sr)


###################################################
### Hazard ratio from weibull model
###################################################

exp(-coef(eg2b.sr)[2]/eg2b.sr$scale)


###################################################
### Plot weibull survival curves
###################################################

curve(pweibull(x,scale=exp(coef(eg2b.sr)[1]),shape=1/eg2b.sr$scale,
               lower.tail=FALSE),from=0, to=max(eg1b$days), ylim=c(0,1),
      ylab="Survival", xlab="Days")
curve(pweibull(x,scale=exp(sum(coef(eg2b.sr))),shape=1/eg2b.sr$scale,
               lower.tail=FALSE),from=0, to=max(eg1b$days), add=TRUE, lty = 2)
lines(survfit(Surv(days,status)~rx,data=eg1b),mark.time=FALSE,
      col="red",lty=c(1,2))


###################################################
### Examining log hazard
###################################################

kmf <- survfit(Surv(days,status)~rx,data=eg1b)
H1 <- -log(kmf[1]$surv)
H2 <- -log(kmf[2]$surv)
Ht1 <- eg1b.km[1]$time
Ht2 <- eg1b.km[2]$time
dH1 <- diff(H1)/diff(Ht1)
dH2 <- diff(H2)/diff(Ht2)
loghaz <- data.frame(ldH=log(c(dH1,dH2)),times=c(Ht1[-1],Ht2[-1]),
                     rx=factor(c(rep(1,length(dH1)),rep(2,length(dH2))),
                       labels=c("Control","Treatment")))

library(lattice)
xyplot(ldH~times,data=loghaz,groups=rx,type=c("p","smooth"),
       ylab="log hazard",xlab="Time")


###################################################
### Fitting log-logistic model
###   To fit log-normal, use dist="lognormal" in
###   the survreg() command below.
###################################################

llog.fit <- survreg(Surv(days,status)~rx,data=eg1b,dist="loglogistic")
summary(llog.fit)

###################################################
### Plotting the result
###################################################

pct <- 0:99/100
llog.pred <- predict(llog.fit,newdata=data.frame(rx=factor(1:2,labels=c("Control","Treatment"))),
                     type="quantile",p=pct)
plot(kmf,mark.time=FALSE,lty=1:2)
matlines(t(llog.pred),cbind(1-pct,1-pct),col="red",lty=1:2)

###################################################
### Adding the Weibul fit for comparison
### This shows an alternative way to obtain the plots
### shown earlier.
###################################################

weib.pred <- predict(eg2b.sr,,newdata=data.frame(rx=factor(1:2,labels=c("Control","Treatment"))),
                     type="quantile",p=pct)
matlines(t(weib.pred),cbind(1-pct,1-pct),col="blue",lty=1:2)